npj Systems Biology and Applications

High-dimensional population balance equations (PBEs) provide a natural framework for modeling heterogeneous cell populations, but their direct numerical solution becomes computationally prohibitive when the internal state space contains many molecular variables. We propose a hybrid mechanistic-machine learning framework for reducing and simulating PBEs defined over high-dimensional intracellular coordinates. The cell population is described by a number density n(x, t), where x [isin] [R]N represents gene and protein states associated with macrophage activation. A dynamics-preserving autoencoder maps this state space to a low-dimensional latent coordinate z [isin] [R]d, with d << N, while retaining key qualitative features of the underlying gene regulatory network, including attractor structure and multistability. Mechanistic information from the original regulatory dynamics is used to construct interpretable drift and diffusion terms for the reduced latent-space PBE. The reduced PBE is solved using a stochastic Lagrangian particle representation, in which particles evolve according to stochastic differential equations (SDEs) corresponding to the latent drift and diffusion fields. The resulting latent-space solution is subsequently decoded and propagated back into the original state space to recover physically interpretable cellular dynamics. We demonstrate the framework on macrophage polarization under cytokine-dependent regulation, including gene knockout perturbations. Overall, the proposed framework provides a computationally tractable and mechanistically interpretable route for integrating single-cell genomic data with population balance models of cell-state dynamics.

Analysis of circadian patterns in time-series data requires computational methods that can accommodate several factors, including variable sampling resolution, replicate number, and missing values. Most existing tools simplify rhythmicity to a strict dichotomy based solely on a single p-value threshold. This leads to a level of uncertainty that affects many biological targets. We developed CODAC (Circadian Oscillation Detection Analysis and Comparison), a framework that integrates nonlinear constrained optimization with a multicriteria rhythmicity classification scheme to evaluate rhythmic patterns without relying on a single statistical cutoff. This approach allows CODAC to identify and exclude medium-confidence rhythms rather than force them into a rhythmic/arrhythmic dichotomy. CODAC comprises four modules: (i) CODAC_single estimates rhythmicity within a single group; (ii) CODAC_flex extends this to identify distinct waveform types within one group; (iii) CODAC_compare performs pairwise comparisons across two or more groups to detect rhythmic or arrhythmic changes; and (iv) CODAC_multi handles more complex designs involving multiple-group comparisons. Using in silico simulations and public transcriptomic datasets, we show that CODAC performs comparably to established methods while providing additional flexibility for rhythm classification and comparison. Taken together, CODAC provides a flexible and open-source package for circadian timeseries analysis with automated visualization tools.

Modern medicine implicitly assumes that physiological responses to intervention are predictably determined by administered treatments. However, physiological systems containing intrinsic delays between the detection of a stimulus and the biological response may violate this assumption. We investigate the human glucose-insulin system as described by the Ultradian model and mathematically demonstrate that clinically relevant forcing protocols-such as pulsatile insulin delivery and step-wise glucose infusion, both commonly used in intensive care units (ICUs)-can induce sustained temporal chaos that may hamper accurate prediction of the physiological response. If not accounted for, these chaotic dynamics could create difficulties in achieving optimal dosing and timing when administering glucose and insulin in clinical or home care settings. This phenomenon, termed delay-induced uncertainty (DIU), arises from the interaction between physiological delay, intrinsic shear near a limit cycle, and external forcing. Using the Ultradian glucose-insulin model, we compute top Lyapunov exponents to quantify predictability. Across a range of pulsatile and step-wise forcing regimes, including stochastic amplitudes drawn from Markov processes, we observe positive Lyapunov exponents, indicating sustained chaos. Our results suggest that delayed endocrine regulation may fundamentally limit the predictive value of the models used to develop glycemic management strategies, with implications for clinical protocols in the ICU.

Cells respond to their environment through protein networks often dysregulated in cancer, making dynamical modelling crucial. Limitations in experimental data and computational resources motivate coarse-graining methods to build low-dimensional descriptions. Yet classical approaches to coarse-grained modelling rely on strong assumptions, leaving it unclear when partial experimental observations support reduced descriptions of system dynamics. Here we show that symbolic regression (SR) provides a data-driven way to test whether, and how compactly, the dynamics of a signalling system coarse-grain over the measured variables, and, when they do, infers mechanistically interpretable models. In synthetic enzyme systems, SR recovers Michaelis-Menten kinetics for the two-step mechanism and under three-step extensions. As data quality is degraded, SR simplifies toward effective kinetic laws while preserving correct theoretical limits. Applied to published time-resolved ERK phosphorylation data, SR identifies compact phospho-ERK rate laws in selected cancer-relevant gene overexpression contexts, yielding interpretable kinetic effects. A sparse neural ODE baseline requires few inputs where SR succeeds, but on average more where it fails, indicating that, where a reduced model is learnable at all, SR failure is associated with more complex dynamics that a simple mathematical model cannot describe. Together, these findings establish symbolic regression as a way to test when a compact coarse-grained description is warranted, generating hypotheses where one holds and motivating potential new measurements where it does not.

Cancer cells rapidly rewire their metabolism, from efficient energy production toward anabolic processes, to sustain uncontrolled growth. Decoding such metabolic shifts is essential for uncovering novel therapeutic targets. To map systems-level metabolic changes across cancer types, we built context-specific genome-scale metabolic models for eight tissues (lung, thyroid, stomach, prostate, liver, kidney, colon, and breast) using gene expression data from The Cancer Genome Atlas (TCGA). Applying constraint-based modelling, we then identified differentially regulated pathways through flux enrichment analysis, revealing tissue-specific rewiring: branched chain amino acid metabolism was suppressed in breast cancer; sphingolipid metabolism was downregulated in colon, kidney, and thyroid but upregulated in breast. We further propose a model-driven pipeline to identify and characterise metabolic vulnerabilities. We first identify synthetic lethal reactions in normal tissues and their corresponding single lethal counterparts in cancers, thereby enabling the identification of metabolic "collateral lethal" reaction pairs for each cancer. Model-predicted collateral lethal gene pairs, including CMPK1-AK in colon, ALDOA-PGD in prostate, and SLC25A26-UQCRB in liver models, were supported through computational validation using DepMap data on gene essentiality. Subsequently, we show how to interpret metabolic rewiring in cancer tissues while accounting for any collateral lethal pairs. In summary, our results establish a systemic framework for decoding metabolic rewiring and synthetic lethal vulnerabilities in cancer.

Neutrophil swarming is a critical immune response in mammals and fish, in which neutrophils are recruited to inflammatory sites where they coordinate into a swarm that neutralizes pathogens. While excessive swarming can drive prolonged inflammation, a quantitative understanding of swarming dynamics remains limited. We developed a one-dimensional radial reaction-diffusion model of neutrophil swarming with two kinetic parameters, in order to capture the self-limiting swarming dynamics in both murine and human neutrophils in response to different inflammatory stimulus sizes. To ensure that the inverse problem is well-posed, we first performed sensitivity and identifiability analyses. We then developed a physics-informed neural network (PINN) to infer the key parameters governing swarm expansion and self-limitation. To account for uncertainty in noisy experimental measurements, we further extended this framework to a Bayesian PINN (B-PINN), which provides credible intervals for the inferred parameters. Both models were validated against synthetic data generated by numerical simulation and subsequently applied to in vitro experimental data from human and murine neutrophils in response to three bioparticle cluster sizes. The PINN-inferred dynamics show that larger bioparticle clusters are associated with greater cumulative recruitment and larger swarms in both species. The models further reveal species-specific differences in both the amplitude of initial recruitment and the timescale on which it self-limits. Additionally, the B-PINN posterior distributions quantify uncertainty in these species- and cluster size-dependent trends and identify where additional measurements would be most informative. To our knowledge, this is the first application of physics-informed machine learning to model neutrophil swarming dynamics. This framework provides a starting point for systematically comparing recruitment dynamics between human and murine neutrophils and offers guidance for future experimental design.

Isotopic Metabolic Flux Analysis (I-MFA) is a standard approach for estimating intracellular metabolic fluxes. I-MFA infers fluxes by comparing simulated and measured metabolite isotopologue distributions (MIDs) of metabolites from isotope labeling experiments. MIDs represent fractional abundances that strictly sum to one for any given metabolite, thus they are inherently compositional data. However, state-of-the-art estimation approaches rely on calculating standard Euclidean distances between MIDs in a non-compositional paradigm, introducing a systemic bias. To resolve this, our study proposes compositional I-MFA. We demonstrate how to construct a meaningful orthonormal basis for MIDs via ordered sequential binary partitioning, which can be used to perform isometric log-ratio (ILR) transformation. As a minimal change to existing I-MFA workflows, we suggest estimating fluxes by minimizing Euclidean distances between ILR-transformed MIDs. We validated this framework against traditional methods using both a toy model and a biologically realistic model, evaluating point estimates, sensitivity across varied true fluxes, and confidence intervals. In the two examples, compositional I-MFA consistently outperformed traditional approaches, reducing mean squared error of flux point estimates by an average of 42.6% and substantially narrowing confidence intervals. We conclude that compositional data analysis significantly improves I-MFA and can be implemented as a simple drop-in replacement for current pipelines. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=156 SRC="FIGDIR/small/742769v1_ufig1.gif" ALT="Figure 1"> View larger version (34K): org.highwire.dtl.DTLVardef@1a53aa4org.highwire.dtl.DTLVardef@ad225aorg.highwire.dtl.DTLVardef@aa430eorg.highwire.dtl.DTLVardef@1880ca_HPS_FORMAT_FIGEXP M_FIG C_FIG HighlightsO_LINew compositional data approach improves metabolic flux estimation. C_LIO_LIThis data transformation requires minimal changes to existing workflows. C_LIO_LIThe new method reduced MSE of flux estimates by 42.6% in two examples tested. C_LIO_LIThe confidence intervals of the estimated fluxes were substantially narrowed. C_LIO_LIEstimation accuracy remained robust across a wide range of metabolic fluxes. C_LI

Mechanistic drug-development programmes often have more biological evidence than they can safely quantify. We developed evidence-constrained mechanistic synthesis (ECMS), a framework that classifies what information each finding contains and converts only that information into restrictions on a family of mechanistic hypotheses. Evidence shifts the frequency of supported events in a reproducible ensemble rather than being converted into unsupported coefficients or probabilities of biological truth. In a chronic spontaneous urticaria (CSU) implementation, a representative, non-exhaustive corpus of 114 atomic findings from 53 sources and 13 public data resources compiled 18 relation/context constraints and a frozen 4,096-hypothesis ensemble. Regimen evaluation was formulated as continuous multi-node target matching: researchers specify desired changes and importance coefficients for modeled nodes, while package-declared controls vary continuously. A deterministic Sobol-to-block-refinement search, validated on all 4,096 hypotheses, reduced target-matching loss by 27.3% relative to the best of 44 deterministic anchors under a prespecified heuristic demonstration profile; changing the objective profile changed the selected control vector without changing the evidence ensemble. A complementary D-only reference analysis localized decision-relevant uncertainty around the mast-cell-to-disease relation, illustrating that mechanistic prioritization depends on the declared objective. ECMS is intended for the pre-calibration stage of drug development: it makes heterogeneous literature computable while keeping evidence, uncertainty and decision preferences distinct.

Long-term in vivo transcriptomic time courses are costly, limiting assessment of chronic molecular responses from short studies. We developed a pharmacokinetics-informed transcriptomic ordinary differential equation model (PKT-ODE) that links an oral pharmacokinetic profile and Hill drug-effect function to first-order turnover of co-expression modules. The model was fitted to rat liver responses to fenofibrate at three doses in Open TG-GATEs through day 8. At the held-out day-29 endpoint, PKT-ODE achieved Pearson r = 0.960 and mean squared error (MSE) = 0.148. In this dataset, these values achieved lower prediction error and higher correlation than four statistical baselines and validation-selected linear and multilayer-perceptron transition models. Literature-curated peroxisome proliferator-activated receptor target genes occurred only in modules with positive fitted drug effects. These results provide a proof of concept for pharmacokinetics-informed transcriptomic extrapolation; cross-compound, cross-organ and alternative-regimen performance remain to be tested.

Crohns disease (CD) is a chronic inflammatory condition of the gastrointestinal tract for which anti-tumour necrosis factor (anti-TNF) agents remain a first-line biologic therapy. However, remission rates are modest, and the mechanistic basis of non-response is poorly characterised. A common resistance mechanism is thought to emerge when alternative inflammatory cascades compensate for TNF inhibition, but the interactions underlying this rewiring have not been systematically characterised. We applied CytokineLink, our previously developed systems immunology framework, to single-cell RNA sequencing data from CD patients sampled before and after anti-TNF therapy. We reconstructed networks of interacting cytokines across samples stratified by treatment phase, response, and inflammation status, and identified condition-specific cytokine interactions and feedback loops, statistically validated against degree-matched random networks. We clustered the generated networks based on their inflammation, response, and treatment status. The pre-treatment inflamed non-responder network contained the largest set of unique interactions, organised around a connected module driven by IL17C targeting downstream TNF, IL6, IL1B, CXCL1/2/3/8, and CCL20. IL17C was produced by a population of non-ileal enteroendocrine cells, differentially abundant at baseline in non-responders. Gene set variation analysis in an independent cohort confirmed elevated non-responder module activity in colonic tissues of non-responders. Feedback loop analysis revealed that responder networks were characterised by persistent IL10 circuits sustained by macrophage populations and acquired tissue-remodelling interactions after therapy, whereas non-responders lost IL10 feedback loops post-treatment and gained TNF-containing motifs, including circuits signalling through the upstream activator TL1A. Our findings characterise the mechanism of anti-TNF non-response as a cytokine network, in which pre-existing epithelial-driven inflammatory modules and the failure to preserve regulatory feedback sustain TNF-independent inflammation in CD. By characterising cytokine interactions at the systems level, our approach moves beyond single-cytokine models of anti-TNF resistance to provide a mechanistic framework for understanding the biological basis of treatment failure in immune mediated diseases.

Control of gene expression by transcription factors (TFs) is a critical mechanism for cells to maintain homeostasis in response to environmental signals. Gene network models that predict regulatory interactions between transcription factors and the genes they control aid in understanding these complex processes. These models are useful as they provide testable hypotheses of regulatory interactions, transcription factor function, and accelerate the study of uncharacterized transcription factors. However, inference of these models is computationally challenging due to the vast quantity of data required given the many possible states of the regulatory network. Microbial genomes encode hundreds of transcription factors, with numerous interactions that require substantial functional genomics datasets to infer. This problem is accentuated in understudied organisms, species that would greatly benefit from an inferred network for biological discovery, where the lack of available data is particularly constraining for effective inference. To address this problem, we have developed GRN-BMuSeR (Gene Regulatory Networks from Bayesian MUlti-SpEcies Regression), a novel multitask approach to gene regulatory network inference that leverages gene orthology between closely related species to improve inference performance. We evaluate its performance on a dataset from the well-studied bacterial species Bacillus subtilis, demonstrating improved performance in multitask settings. Applying the model to simulated data reveals utility in multi-species contexts. Finally, we apply our models to infer GRNs and explore predictions for two hypersaline-adapted archaeal species. We leverage a rich dataset from Halobacterium salinarum to inform the inference of the gene regulatory network of Haloferax volcanii, for which a more limited genomics dataset was available. We generate a large compendium of gene expression data for Hfx.volcanii for GRN inference input. Through exploration of resultant network predictions, we show concordance with known TF functions and discover hundreds of novel TF functional predictions. Moving forward, our results provide a framework to generate testable hypotheses that will serve to guide experimental work and accelerate discovery in these understudied species.

Microorganisms interact through the exchange of metabolites and competition for shared substrates, and this metabolic coupling shapes the composition and function of microbial communities. Community flux balance analysis (cFBA) can predict such behaviour - the maximum community growth rate, the metabolic fluxes and the relative abundances of the species - from stoichiometric models of their metabolism, but existing formulations are either complex and hard to scale as communities grow or cannot predict optimal growth rates. Here we present a physiology-based formulation of cFBA in which each species' metabolism is reduced to a few macrochemical equations, one for each 'metabolic mode' the species can use, and the whole community is then solved as a single linear program. From this, the method predicts the optimal composition of the community, its maximum growth rate, the metabolites exchanged between the species, and the net conversion the community carries out as a whole; its ecological service. This reduction makes it far simpler to build and solve models of larger communities. We illustrate the approach on a two-species synergistic community that can be verified by hand, apply it to a five-member anaerobic digestion community, and use it to predict the metabolic interactions of a genome-scale syngas-fermenting coculture. Characterising these communities at their optimal steady states, we show that each species is driven to a distinct metabolic strategy. We discuss the method both as a practical tool for larger microbial communities and as a means of uncovering the ecological principles that govern them.

Light simultaneously provides phototrophic organisms with energy and with information about environmental time. These two functions need not impose the same response to fluctuations in irradiance: photosynthetic outputs should remain amplitude-sensitive, whereas circadian phase should reject changes that do not alter dawn, dusk, or photoperiod. We formulate this problem for two spectral inputs by decomposing their logarithmic intensities into a common-irradiance coordinate a and a spectral-contrast coordinate r. The contribution of channel i to phase is Qi = ZiGi, where the non-negative gate Gi determines when the pathway is active and the signed phase-response projection Zi determines whether this activity advances or delays the oscillator. For a locked oscillator, robustness to common irradiance together with retained contrast sensitivity requires two non-zero cycle-averaged contributions of opposite sign, A1 [~=] -A2 = 0. Energetic responsiveness is preserved only when the physiological projection of the same inputs is not proportional to their phase projection. A canonical repressilator provides an explicit nonlinear realization of these conditions. Positive gates placed on opposite lobes of its infinitesimal phase-response curve strongly attenuate common-mode phase shifts while preserving contrast sensitivity. A minimal photosynthetic-capacity model then shows how this organization protects temporal alignment under day-to-day irradiance fluctuations. At the largest variability tested, differential routing reduced the mean phase displacement by more than one half and the associated alignment loss by approximately 82%, whereas the resulting production advantage remained small, approximately 0.1%. Thus, multichannel light sensing can stabilize circadian timing without suppressing the energetic response to irradiance. HighlightsO_LIAnalytical routing conditions separate common irradiance from spectral contrast. C_LIO_LIPositive temporal gates can generate opposite signed phase contributions. C_LIO_LIPhase robustness requires a projection distinct from the energetic projection. C_LIO_LIA canonical oscillator provides a constructive illustration of the mechanism. C_LIO_LIThe functional benefit is improved temporal alignment rather than a large growth gain. C_LI

Gene regulatory network models treat interaction parameters as fixed, although regulatory efficacy fluctuates. We asked how temporal fluctuations in interaction strength reshape phenotype occupancy in cell-fate decision GRN motifs. Across large parameter ensembles, anchored fluctuations largely preserved deterministic occupancies. Additive fluctuations increased occupancy of all-high co-expression states, particularly where high expression was accessible. In contrast, multiplicative fluctuations biased inhibitory interactions toward stronger repression and favored single-high states in a topology-dependent manner. Deterministic controls sampled from noise-induced parameter distributions did not fully reproduce these effects. A Boolean-limit analysis revealed an intrinsic upward bias: loss of repression increased expression regardless of regulator state, whereas stronger repression acted only when the regulator was present. Analyses of epithelial-mesenchymal plasticity and gonadal-fate networks showed increased occupancy of hybrid team-expression states under additive fluctuations. Thus, regulatory noise can reshape the developmental landscape in opposing directions, pushing cell-fate systems toward either progenitor-like or terminally differentiated states.

Cells coordinate diverse biological processes through interactions among thousands of molecules, but mapping these interactions comprehensively and systematically remains an open problem. Pathway reconstruction algorithms address this problem by linking molecules of interest, identified from high-throughput omics experiments, using prior knowledge encoded as background interaction networks. This process recovers intermediate molecules and interactions that were not directly measured in the experimental data but plausibly connect the observed molecules. It generates testable hypotheses about interactions that drive cell behavior and informs the choice of follow-up experiments. Many algorithms have been created over decades, each optimizing different computational objectives and relying on different assumptions. The resulting heterogeneity has made benchmarking challenging, limiting systematic comparisons. Therefore, selecting an algorithm for a given biological context remains a non-trivial and poorly informed task. This registered report presents a large-scale benchmark of pathway reconstruction algorithms, evaluating 14 algorithms across 822 datasets from four biological settings. To enable this benchmark, we introduce Signaling Pathway Reconstruction Analysis Streamliner (SPRAS), which standardizes algorithm inputs, outputs, and execution into a formal framework, enabling systematic comparison that was previously infeasible. We will assess each algorithm on reconstruction performance against gold standard pathways, algorithm similarity, and computational performance across different biological contexts. Together, these evaluations will provide quantitative evidence for understanding pathway reconstruction algorithm behavior and guiding algorithm selection.

Background and Objective: Predicting treatment outcomes in infectious diseases requires accounting for the interplay between drug effects, pathogen dynamics, and host immunity. Integrating pharmacological and immunological approaches into a single simulation environment remains a fundamental challenge in both theory and practice. We aimed to develop and validate a multiscale in silico framework coupling these processes, and to quantify their respective contributions to bacterial clearance. Methods: We present the Drug-Host-Pathogen Interaction (DHPI) framework, combining three independent mechanistic components: a physiologically based pharmacokinetic model of drug disposition, a pharmacokinetic-pharmacodynamic model of drug-induced bacterial killing, and a stochastic agent-based model of the immune response. Continuous concentration profiles are time-averaged onto the agent-based time grid, assigned to bacterial phenotypic states, and converted into per-agent killing probabilities, so that drug-mediated and immune-mediated death events are recorded separately at each step. The framework was applied to simulate symptomatic pulmonary tuberculosis. Phenotype-specific drug-efficacy parameters were inferred using Approximate Bayesian Computation from historical clinical data on eight weeks of 600 mg rifampicin monotherapy, and validated against independent early bactericidal activity data over a disjoint time window. Results: The calibrated framework reproduced the observed decline in bacterial load, and matched reported early bactericidal activity over the first week. In a virtual cohort of symptomatic patients, drug-mediated killing accounted for 81-88% and immune-mediated killing for 12-19% of total bacterial elimination over the 60-day treatment course, while the dormant, granuloma-contained fraction rose from 0.20-0.29 in the first week to 0.85-0.89 at treatment completion. Over a follow-up of up to 50 years, patients reaching clinical cure had accumulated more memory lymphocytes during treatment than those progressing to clinical failure or death; moreover, the final outcome depended on the immune changes occurring during therapy rather than on the initial disease stage. Conclusions: The results show that the DHPI framework can reproduce treatment dynamics observed in patients and enable the analysis of how therapy reshapes host immune responses and subsequent disease trajectories. By explicitly representing drug-host-pathogen interactions, it provides a mechanistic basis for in silico treatment simulations and for the study of long-term immune consequences of antimicrobial therapy.

Tissue shape and function are defined by the mechanical interactions of cellular and extracellular components. Lymph nodes cyclically remodel in response to immune challenges whilst preserving essential stromal structures. However, the relative contributions of the fibroblastic reticular stromal cell network and the ensheathed extracellular matrix, remain undefined. We quantified the contribution of ECM to the viscoelastic properties of lymph nodes to parameterise an in silico model exploring the FRC network's adaptation to pressure-driven tissue expansion. The balance between tissue pressure, FRC contractility and ECM stiffness permit robust remodelling and growth, while maintaining physiological geometries and balancing force distribution. Local perturbation of ECM stiffness or FRC contractility disrupts force distribution globally and impacts FRC proliferation and tissue expansion. Spatially dispersed perturbations exert higher impact on tissue architecture than equivalent localised perturbations, with effects propagating across the network. The integration of cellular and extracellular mechanics thereby enables robust lymph node remodelling.

Zero-dimensional chemical master equations, ordinary differential equations, and compartmental population models replace spatial stochastic biological systems by vectors of total counts or densities. This study asks when that projection is exact and whether information retained in spatial correlations can diagnose its practical failure. Exact Markov closure is characterized by an aggregate-rate lumpability condition: for every retained transition, the sum of microscopic transition rates must be constant over all spatial configurations with the same counts. Violations are connected to BBGKY-type correlation hierarchies and to mean-field, pair, and triplet closures. Conditional rate, finite-time predictive, memory, path-space, and correlation Kullback-Leibler risks quantify distinct losses. An exactly solvable two-compartment reaction separates structural non-closure from recovery of a well-mixed law under fast hidden mixing. Copy number and a spatial mixing-interaction ratio connect concentration, volume, diffusion, and reaction parameters to practical screening, including an Escherichia coli-scale example. The same projection logic is evaluated in controlled spatial susceptible-infectious-removed and predator-prey benchmarks. Across mixed and segregated initial conditions and four mobility regimes, pair-correlation risk was strongly associated with the error of the corresponding zero-dimensional ordinary differential equations (Spearman correlations 0.95 and 1.00; pooled 0.99). A nearest-neighbour exchange sensitivity analysis preserved the positive risk-error ranking. These benchmarks do not establish a universal threshold, but support correlation information as a transferable diagnostic for selecting among count, pair, higher-order, and explicit spatial descriptions.

Although most human protein coding genes have functional annotations in databases, such as GeneCards, many remain poorly characterized. To address this gap, computational tools can be leveraged to predict the functional roles of under-annotated genes by extracting patterns from complex biological networks. Here we introduce Brain-for-Biotech (BfBio), a framework designed to identify genes important for vascular endothelial cells (EC), which are crucial cells for vessel formation (angiogenesis), vascular homeostasis, hemostasis and blood/tissue barrier function but also critical mediators of immunity and cancer progression. BfBio utilizes a Personalized PageRank (PPR) algorithm on an integrated network of different omics datasets and publicly available gene-gene/protein-protein interaction databases. In this study, we apply the predictive capabilities of BfBio to infer angiogenic stalk cell phenotype function in genes for which this function was not known before. By leveraging a set of genes characterizing the stalk cell cluster in lung tumor EC models previously identified, we have achieved a high Area Under Receiver Operative Characteristic (AUC-ROC) performance (0.837). Enrichment analysis, coupled with a text mining application, further confirmed that among the 49 predicted genes four of them were poorly characterized yet possessed biologically relevant properties and were linked to cancer, thereby validating BfBio as a robust tool for prioritizing novel therapeutic targets in vascular biology.

We present BIOCURRENT, a causal inference framework that reconstructs donor-specific pseudotime geometry in transcriptomic data. By modeling gene expression as a function of baseline characteristics, microenvironmental context, and latent pseudotime, BIOCURRENT enables comparison of compressed or expanded progression intervals across transcriptional state transitions. We introduce $\Delta\Delta T$, a geometry-based estimator that quantifies differences in pseudotime intervals across conditions, enabling evaluation of changes in pseudotime intervals under hypothetical modulation of microenvironmental programs. Applications to thymic T-cell developmental lineages and to COVID-19 immune dysregulation reveal condition- and donor-specific distortions of progression intervals. Counterfactual simulation links microenvironmental context to changes in specific intracellular state transition intervals. By localizing deviations in pseudotime geometry, BIOCURRENT identifies whether shifts in transcriptomic programs emerge early or later along transcriptomic coordinates and reveals upstream programs associated with these distortions. Such localization supports transcriptional stage-aware mechanistic hypotheses and suggests candidate intervention checkpoints in complex biological systems.